www.gusucode.com > 基于Matlab的MIMO通信系统仿真 含报告；司中威；了解移动通信 > 基于Matlab的MIMO通信系统仿真 含报告；司中威；了解移动通信关键技术，了解数字通信系统仿真流程，实现基本的信道编译码、调制解调等通信模块。（好评如潮，课设拿满） 学习并实现MIMO空时处理技术 学习性能分析的思路和方法/mimo/matlab for mimo 2x2/rx.m

     function BER = rx(detect_type,select_type)
% Transmission
% Maxime Maury
% 05-04-25

% If argument is omitted, the channel is simulated
if (nargin < 1)
    detect_type =3;
end
if (nargin < 2)
    select_type = 1;
end

close all;
disp('-- Start RX --')

% ---------------------------------------------------
% Global parameters
% ---------------------------------------------------

disp('Loading...');

% Real or simulated channel?
load('SimType');

% Load global parameters
load('tx_param');

% Time bewteen 2 switches within a frame (before upsampling)
switch_len = guard_len + training_s_len + training_len;

% Delta Symbols set for exhaustive searching
[Delta_S_Set, Set_len] = D_Symbol_Set;

% Initial Antenna selection from {1,2,3,4} but only 4 different choices:
% that are: 1&3, 1&4, 2&3, 2&4
A = 1;
B = 3;

% Type of different detectors, parameters for Detector.m 
ML = 1;
JMMSE = 2;
ZF = 3;

% Type of different antenna selection criteria methods 
MBER = 1;
MMI = 2;
LAZY = 3;
MNP = 4; %Minimum Noise Power
MMNP = 5;
LAZY2 = 6;
% (De) activate antenna selection
antenna_select = 1;
if (real_ch==1)
    antenna_select = 0;% disable antenna selection when use real channel
end


% H channel matrices Initial
H = zeros(4,nr_frames*2+2);
H([A B],1:2) = eye(2); %Initial
H([3-A 7-B],1:2) = eye(2); 

% H channel matrices Initial (simulated channel)
if (real_ch==0)
    H_s = zeros(2,frames_len*2);
end

% Memory for the channel estimation
lambda_H = 0.95; 

% Load the channel applied 
if (real_ch==0)
    load('ChannelMatrix')
end

% Read data file (binary data sent via the TX)
fid = fopen('transmit.dat', 'r');
total_data_len = data_len_bits*2*nr_frames;
data_sent = fread(fid, total_data_len, 'int16');
data_sent = data_sent';
fclose(fid);

% Read output of the channel
load('ChannelOutput');

% data_rec is the data received in bits
data_rec = zeros(1,data_len_bits*2*nr_frames);

% data_rec_split is the data received in bits on Channel 1 and on Channel 2
data_rec_split = zeros(2,data_len_bits*nr_frames);

% M_sync is the analysis window length for synchronization
M_sync = 250;
    
% Out_match is the output of the match filter
Out_match = zeros(4,Lf+M_sync);

% Out_data data read in one frame
Out_data = zeros(4,data_len); % Data in symbols

% For plotting the constellation
r1_data_I_block = [];
r1_data_Q_block = [];

% ---------------------------------------------------
% Noise estimation
% ---------------------------------------------------

disp('Nois estimation...');

% Estimate the noise variance from the first values
sigma_sqr_noise = (var(Y_R(A,1:init_len)) + var(Y_R(B,1:init_len)) )/2;
sigma_sqr_noise_bb = 4 * sigma_sqr_noise;
N0 = sigma_sqr_noise_bb/2/L;
Var_symbol =  (var(Y_R(A,2000:2000+init_len)) + var(Y_R(B,2000:2000+init_len)) )/2;
Var_symbol = Var_symbol - sigma_sqr_noise;
Es_bb = Var_symbol * norm(pulse_shape)^2/L;
SNRbb = 10 * log10(Es_bb/sigma_sqr_noise_bb);
SNRVar = 10 * log10(Var_symbol/sigma_sqr_noise);

disp(['Estimated bb noise variance: ', num2str(sigma_sqr_noise_bb)]);
disp(['Var(Signal)/Var(Noise):      ', num2str(SNRVar)]);

% For plots
Nc_sync = ceil(sqrt(nr_sync_frames));
Nr_sync = ceil(nr_sync_frames/Nc_sync);
Nc_est = ceil(sqrt(nr_frames*2));
Nr_est = ceil(nr_frames*2/Nc_est);
Nc = 3;
Nr = 2;

% Length of the signal
signal_len = size(Y_R,2);

% ---------------------------------------------------
% Start of the signal detector
% ---------------------------------------------------

disp('Detection of the signal...');

% Determine the beginning of the signal
sig_start = 0;

% Magnitude of the signal at the beginning
signal_mag = sigma_sqr_noise_bb;

% Magnitude threshold of signal detected
mag_T = 1.2*sqrt(Var_symbol); 
if SNRVar > 10
    mag_T = 1.0*sqrt(Var_symbol);
end
% Pourcent of the last value kept
lambda = .95; 

values = [];
while ( signal_mag < mag_T )

    sig_start = sig_start + 1;
    
    % If we are at the end the detection has failed
    if (sig_start>signal_len)
        error('Signal detection has failed')
        break
    else
        signal_mag = lambda*signal_mag + (1-lambda) * abs(Y_R(A,sig_start));
        values = [values signal_mag];
    end
end

disp(['The signal starts at instant:',num2str(init_len+1)]);
disp(['Estimation:                  ',num2str(sig_start)]);

% Excepted time instant when the 1st frame starts
t_frame1 = sig_start + time_end + 1;

% The signal detection algorithm will not detect the very begining of
% the sinusoid. Instead, a delay will be genereated. Then, we need to
% compensate for this delay and stop the DPPL in advance:
estimated_delay = 100;
DPLL_end = time_end - estimated_delay;

% Plot the estimated times
% Plot the beginning of the signal
figure(1);
subplot(Nr,Nc,1)
plot(Y_R(A,1:init_len + time_end + 500)); hold on;
plot(sig_start,0,'pg'); hold on;
plot(DPLL_end + sig_start, 0 ,'sy'); hold on;
plot(t_frame1,0,'*k');
legend('Signal','Detected start','End of DPLL','Frame')

% Plot the signal detection output
subplot(Nr,Nc,2);
plot(values)
hold on;
plot(mag_T*ones(1,length(values)),'r');
title('Signal detection');
grid on;

% ---------------------------------------------------
% DPLL
% ---------------------------------------------------

subplot(Nr,Nc,3);
hold on;
[deltaA, thetaA] = lock_phase(2*pi*(Fc/Fs), 0, 0, Y_R(A,sig_start:DPLL_end+sig_start));
[deltaB, thetaB] = lock_phase(2*pi*(Fc/Fs), 0, 0, Y_R(B,sig_start:DPLL_end+sig_start));
title('Frequency offset estimation: Output of the DPLL');
grid on;
drawnow

% %test the zero-crossing method
% Fc_1 = frequency_estimation(Y_R(A,sig_start:DPLL_end+sig_start),Fs);
% Fc_2 = frequency_estimation(Y_R(B,sig_start:DPLL_end+sig_start),Fs);
% Fc_zero = (Fc_1 + Fc_2)/2

freq_offset_hat = (deltaA + deltaB)/(2*2*pi);
disp(['Estimated frequency:         ', num2str(Fc+freq_offset_hat*Fs)]);
Fc_hat = Fc + freq_offset_hat*Fs;

for fr=1:nr_frames
    
    n_fr_sync = floor((fr-2)/refresh_fr) + 1;
    n_fr_reg =  (fr-1) - n_fr_sync;
    
%     disp(' ');
%     disp(['Frame ',num2str(fr),' Using antennas ', num2str(A),num2str(B)]);
    
    % Take a bunch of values selected by a window
    % All the following times are used for windowing. They are not expected
    % to be accurate
    t_start =   t_frame1 + n_fr_reg * Lf + n_fr_sync * (Lf+training_s_len*L) ; % Position in the frame  (in symbols)
    t_switch =  t_start + switch_len*L + M_sync;
    
    % Switch to the complementary set of the best set of antenna
    % Worst set selection (I=3-A, J=7-B)
    
    % Demodulation: Downconversion and Pulse-shaping
    
    if (antenna_select)
        I = 3 - A;
        J = 7 - B;
    else
        I = A;
        J = B;        
    end
    
    % ---------------------------------------------------
    % Demodulation and pulse-shaping
    % ---------------------------------------------------

    [rI_I, rI_Q, rJ_I, rJ_Q] = demodulate(Y_R(:,t_start:t_switch),I,J,Fc_hat, Fs, t_start, pulse_shape);        

    % ---------------------------------------------------
    % Synchronize
    % ---------------------------------------------------
    
    if (mod((fr-1),refresh_fr) == 0)
        
        % First synchronization 
%         figure(2);
%         subplot(Nr_sync,Nc_sync, floor((fr-1)/refresh_fr) + 1)
        t_first_train = synchronize(    rI_I, rI_Q, ...
            rJ_I, rJ_Q, ...
            tr_sync1_I, tr_sync1_Q , ...
            M_sync, L , 0  );
        
%         title(['Synchronization - Frame ', num2str(fr)]);
%         drawnow
        
        % Where the first sychronization training sequence stops 
        t_first_train_e  = t_first_train + training_s_len * L - 1;     
    else
        t_first_train_e  = t_first_train - 1;  
    end
          
    
    % Switch at t_first_train_e to the best set! (A, B)
    
    % ---------------------------------------------------
    % Channel estimation
    % ---------------------------------------------------
 
    % Start of the channel estimation training sequence 
    i_sta = t_first_train_e + 1;
    
     % End of it
    i_end = i_sta + training_len*L - 1 ;
 
%     figure(3);                           
%     subplot(Nr_est,Nc_est,2*fr-1)
    
    H_hat_1 = channel_estimation(   rI_I(i_sta:L:i_end), ...
                                    rI_Q(i_sta:L:i_end), ...
                                    rJ_I(i_sta:L:i_end), ...
                                    rJ_Q(i_sta:L:i_end), ...
                                    tr1_I, tr1_Q , tr2_I, tr2_Q ...
                                 );
                                
%     title(['Ch. estimation 1 - Frame ', num2str(fr)]);
%     drawnow

  
    % Best set selection
      
    % All the following times are used for windowing. They are not expected
    % to be accurate
    
    if (mod((fr-1),refresh_fr) == 0)
        t_start = t_start + switch_len*L;
    else
        t_start = t_start + (switch_len-training_s_len)*L;  
    end
    
    t_end = t_start + (data_len + guard_len + training_len)*L + M_sync;
    
    % Demodulation
    [rA_I, rA_Q, rB_I, rB_Q] = demodulate(Y_R(:,t_start:t_end),A,B,Fc_hat, Fs, t_start, pulse_shape);
                        
    % Channel estimation
    
    % Start of the channel estimation training sequence 
    i_sta = t_first_train ;
    
    % End of it
    i_end = i_sta + training_len*L - 1;

%     figure(3);
%     subplot(Nr_est,Nc_est,2*fr)  
    
    H_hat_2 =  channel_estimation(  rA_I(i_sta:L:i_end), ...
                                    rA_Q(i_sta:L:i_end), ...
                                    rB_I(i_sta:L:i_end), ...
                                    rB_Q(i_sta:L:i_end), ...
                                    tr1_I, tr1_Q , tr2_I, tr2_Q ...
                                 );

                                 
%     title(['Ch. estimation 2 - Frame ', num2str(fr)]);
%     drawnow
                                
    
    % Update H estimations
    
    % Previous estimation Hr_n-1
    H_pre_1 = H([3-A 7-B],fr*2-1:fr*2);
     
    % Update the estimation by Hr_n = lambda*Hr_n-1 + (1-lambda)*H_n
    if fr > 1
        H([3-A 7-B],fr*2+1:fr*2+2) = lambda_H*H_pre_1 + (1-lambda_H)*H_hat_1;
    else
        H([3-A 7-B],fr*2+1:fr*2+2) = H_hat_1;
    end
 
    % Previous estimation Hr_n-1
    H_pre_2 = H([A B],fr*2-1:fr*2);
    
    % Update the estimation by Hr_n = lambda*Hr_n-1 + (1-lambda)*H_n
    if fr > 1
        H([A B],fr*2+1:fr*2+2) = lambda*H_pre_2 + (1-lambda)*H_hat_2;
    else
        H([A B],fr*2+1:fr*2+2) = H_hat_2;
    end   
    
    if (real_ch==0)
        
        if (mod((fr-1),refresh_fr) == 0)
            actual_len = Lf_sync;
        else
            actual_len = Lf;
        end
        
        posi = 2*(n_fr_reg * Lf + n_fr_sync* (Lf+training_s_len*L));
        
        for k=1:actual_len
            H_s(:, posi+2*k-1:posi+2*k) = H_hat_2;
        end    
    end
    
    % Downsample
    % Start of the data in the frame
    % Use synchronization from the second sychronization training sequence
    data_start = i_end + 1;
    data_end = data_start + L*data_len - 1;
    
    rA_data_I = rA_I(data_start:L:data_end);
    rA_data_Q = rA_Q(data_start:L:data_end);
    rB_data_I = rB_I(data_start:L:data_end);
    rB_data_Q = rB_Q(data_start:L:data_end);
   
    %%%%%% end of Channel estimation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
          
        
    % ---------------------------------------------------
    %  Antenna Subset Selection, by type = MBER,MMI,MNP,MMNP
    % ---------------------------------------------------
 
    if (antenna_select)
        [A, B]= antenna_selection(H(:,fr*2+1:fr*2+2),Delta_S_Set,select_type,SNRVar);
%          disp(['Frame ', num2str(fr),'Antenna Selection=',num2str(A),num2str(B)])

%         %-----------------verify the DSP code with MATLAB--------------------------        
%         write_float32_hex('C:\ti\myprojects\mimo\receiver\Hr.dat',real(H(:,fr*2+1:fr*2+2))');
%         write_float32_hex('C:\ti\myprojects\mimo\receiver\Hi.dat',imag(H(:,fr*2+1:fr*2+2))');
%         write_float32_hex('C:\ti\myprojects\mimo\receiver\SNR.dat',10^(SNRVar/10));
%         SNR=read_float32_hex('C:\ti\myprojects\mimo\receiver\SNR.dat');

    end
        
    % Detector using type = ML,JMMSE, or ZF to find S_hat
    [rA_data_I, rA_data_Q, rB_data_I, rB_data_Q] =...
        Detector(H_hat_2,rA_data_I, rA_data_Q, rB_data_I, rB_data_Q,detect_type,SNRVar);

    %%%%%%%%%%%%%%%%Minimum distance de-modulator%%%%%%%%%%%%%%%%%%%%%%
    
    % For plotting the constellation
    r1_data_I_block = [r1_data_I_block, rA_data_I];
    r1_data_Q_block = [r1_data_Q_block, rA_data_Q];
    
    % Start of the data burst
    s_start = 1 + data_len_bits*(fr-1);
    
    % End of the data burst
    s_end = data_len_bits*(fr);

    % Detector from symbols to bits
    data_split_rec(1,s_start:s_end) = detect(rA_data_I , rA_data_Q , 1 );
    data_split_rec(2,s_start:s_end) = detect(rB_data_I , rB_data_Q , 1 );

end

% Parallel To Serial
data_rec(1:2:end) = (data_split_rec(1,:));
data_rec(2:2:end) = (data_split_rec(2,:));

% Find where the errors are
n_bits = data_len_bits*2*nr_frames;
positions1 = zeros(1,data_len_bits);
positions2 = zeros(1,data_len_bits);
for p=1:n_bits
    if(data_rec(p) ~= data_sent(p))
        p = mod(p, data_len_bits*2)+1;
        if mod(p,2)
            positions1((p+1)/2)=positions1((p+1)/2)+1;
        else
            positions2(p/2)=positions2(p/2)+1;
        end
    end
end

% Plot the location of the errors in a frame
figure(1);
subplot(Nr,Nc,4);
stem([1:data_len_bits],positions1)
title('Bit errors on stream 1')
subplot(Nr,Nc,5);
stem([1:data_len_bits],positions2,'r')
title('Bit errors on stream 1')

% Compute the BER

BER = sum(abs(data_rec-data_sent))/(length(data_sent));
disp(['BER = ', num2str(BER)]);
disp('End');

% ---------------------------------------------------
% Training sequence constellation
% ---------------------------------------------------

figure(1);
subplot(Nr,Nc,6);
plot(rA_I(i_sta:L:i_end),rA_Q(i_sta:L:i_end),'*r');
grid on;
hold on;
plot(tr1_I,tr1_Q,'*b')
axis equal
title('Constellation of the ch. training seq')

% ---------------------------------------------------
% Channel estimation results
% ---------------------------------------------------

figure(4);
Nc_c = 2;
Nr_c = 1;

if (real_ch)
    
    subplot(Nr_c,Nc_c,1);
    hold on;
    grid on;
    plot(angle(H(1,1:2:end)),'b');
    plot(angle(H(1,2:2:end)),'b-. ');
    plot(angle(H(2,1:2:end)),'g');
    plot(angle(H(2,2:2:end)),'g--');
    plot(angle(H(3,1:2:end)),'r');
    plot(angle(H(3,2:2:end)),'r-. ');
    plot(angle(H(4,1:2:end)),'k');
    plot(angle(H(4,2:2:end)),'k--');
    legend('h11','h12','h21','h22','h31','h32','h41','h42')
    ylabel('Phase');
    
    subplot(Nr_c,Nc_c,2);
    hold on;
    grid on;
    plot(abs(H(1,1:2:end)),'b');
    plot(abs(H(1,2:2:end)),'b-. ');
    plot(abs(H(2,1:2:end)),'g');
    plot(abs(H(2,2:2:end)),'g--');
    plot(abs(H(3,1:2:end)),'r');
    plot(abs(H(3,2:2:end)),'r-. ');
    plot(abs(H(4,1:2:end)),'k');
    plot(abs(H(4,2:2:end)),'k--');
    legend('h11','h12','h21','h22','h31','h32','h41','h42')
    ylabel('Amplitude');
    
else
    subplot(Nr_c,Nc_c,1);
    hold on;
    grid on;
    plot(angle(H_s(1,1:2:end)),'b');
    plot(angle(H_s(1,2:2:end)),'b-. ');

    plot(angle(H_s(1,1:2:end)),'r');
    plot(angle(H_s(2,2:2:end)),'r-. ');

    legend('h11','h12','h31','h32')
    ylabel('Phase');
    title('Estimates')
    
    subplot(Nr_c,Nc_c,2);
    hold on;
    grid on;
    
    plot(abs(H_s(1,1:2:end)*2),'b');
    plot(abs(H_s(1,2:2:end)*2),'b-. ');
    plot(abs(H_s(2,1:2:end)*2),'r');
    plot(abs(H_s(2,2:2:end)*2),'r-. ');
    
    plot(abs(H_var(1,1:2:end)),'m');
    plot(abs(H_var(1,2:2:end)),'m-.');
    plot(abs(H_var(3,1:2:end)),'g');
    plot(abs(H_var(3,2:2:end)),'g-.');    
    
    legend('h11','h12','h31','h32','h11','h12','h31','h32')
    ylabel('Amplitude');    
    title('Estimates and Real values ')
end

% ---------------------------------------------------
% Constellation
% ---------------------------------------------------

if (detect_type~=ML)
    figure;
    Nr = 1;
    Nc = 1;
    
    subplot(Nr,Nc,1);
    plot(r1_data_I_block,r1_data_Q_block,'*b')
    grid on;
    title('Constellation')
end